% computation of stress and force fields corresponding to
% all given functions needed for an exact solution
% of an elasticity problem in Laplace domain:
% -div sigma + rho s^2 u = f
%
% Last Modified: November 21, 2017


div=@(x,y,z) ux(x,y,z)+vy(x,y,z)+wz(x,y,z);
divx=@(x,y,z) uxx(x,y,z)+vxy(x,y,z)+ wxz(x,y,z);
divy=@(x,y,z) uxy(x,y,z)+vyy(x,y,z)+wyz(x,y,z);
divz=@(x,y,z) uxz(x,y,z)+vyz(x,y,z)+wzz(x,y,z);

sxx=@(x,y,z) 2*mu(x,y,z).*ux(x,y,z)+lam(x,y,z).*div(x,y,z);
syy=@(x,y,z) 2*mu(x,y,z).*vy(x,y,z)+lam(x,y,z).*div(x,y,z);
szz=@(x,y,z) 2*mu(x,y,z).*wz(x,y,z)+lam(x,y,z).*div(x,y,z);

sxy=@(x,y,z) mu(x,y,z).*(vx(x,y,z)+uy(x,y,z));
sxz=@(x,y,z) mu(x,y,z).*(wx(x,y,z)+uz(x,y,z));
syz=@(x,y,z) mu(x,y,z).*(wy(x,y,z)+vz(x,y,z));

fx=@(x,y,z) -(lam(x,y,z).*divx(x,y,z)+lamx(x,y,z).*div(x,y,z)...
    +2*mu(x,y,z).*uxx(x,y,z)+2*mux(x,y,z).*ux(x,y,z)...
    +mu(x,y,z).*(vxy(x,y,z)+uyy(x,y,z))+muy(x,y,z).*(vx(x,y,z)+uy(x,y,z))...
    +mu(x,y,z).*(wxz(x,y,z)+uzz(x,y,z))+muz(x,y,z).*(wx(x,y,z)+uz(x,y,z)))...
    +rho(x,y,z).*s.^2.*u(x,y,z);

fy=@(x,y,z) -(lam(x,y,z).*divy(x,y,z)+lamy(x,y,z).*div(x,y,z)...
    +2*mu(x,y,z).*vyy(x,y,z)+2*muy(x,y,z).*vy(x,y,z)...
    +mu(x,y,z).*(wyz(x,y,z)+vzz(x,y,z))+muz(x,y,z).*(wy(x,y,z)+vz(x,y,z))...
    +mu(x,y,z).*(uxy(x,y,z)+vxx(x,y,z))+mux(x,y,z).*(uy(x,y,z)+vx(x,y,z)))...
    +rho(x,y,z).*s.^2.*v(x,y,z);

fz=@(x,y,z) -(lam(x,y,z).*divz(x,y,z)+lamz(x,y,z).*div(x,y,z)...
    +2*mu(x,y,z).*wzz(x,y,z)+2*muz(x,y,z).*wz(x,y,z)...
    +mu(x,y,z).*(uxz(x,y,z)+wxx(x,y,z))+mux(x,y,z).*(uz(x,y,z)+wx(x,y,z))...
    +mu(x,y,z).*(vyz(x,y,z)+wyy(x,y,z))+muy(x,y,z).*(vz(x,y,z)+wy(x,y,z)))...
    +rho(x,y,z).*s.^2.*w(x,y,z);